With the development of ultrafast optics, THz wave has become a powerful tool in studying the nonlinear optical responses of materials at THz frequency^[1–6] for acquiring extensive information about the carrier and quasi-particle dynamics in semiconductors. For example, the coherent THz emission of n-type GaAs,^[1] carrier-wave Rabi oscillations^[2] and coherent ballistic transport of electrons^[3] have been excited using intense THz radiation. There has been observed nonlinear THz absorption bleaching due to intervalley scatterings in Ge, Si, and GaAs.^[4] Razzari et al. studied nonlinear ultrafast modulations of the optical absorption of intense few-cycle terahertz pulses in n-doped InGaAs films recently.^[5] The studies of the nonlinear optical response and dynamic mechanisms have important theoretical meanings, and can provide new data to improve the performances of materials and devices.

InSb and InAs semiconductors are becoming important in many optoelectronic and electronic device applications^[7–10] because of their high mobilities and narrow band gaps. The two materials are commonly used as a traditional THz wave source when irradiated with femtosecond laser pulses.^[11–14] There are also lots of researches on the nonlinear optical response of narrow band gap material in the THz fields. Scholars often use the samples with fixed doping levels excited by the intense THz fields without any direct optical electron–hole generation. For instance, the nonlinear absorption is observed in InSb, and the impact ionization process has been resolved by THz-pump/THz-probe (TPTP) measurements on a femtosecond time scale;^[15,16] the reflection geometry is used to study the intervalley scattering and impact ionization in InAs with intense terahertz pulses.^[17]

In spite of the numerous nonlinear studies already performed, the nonlinear THz radiation under the intense THz field has not been extensively investigated. Moreover, in most of these papers the classical drift diffusion model and Drude model were used for theoretical analysis, and the extremely complex transient process of carriers is difficult to solve accurately. The ensemble Monte Carlo simulation^[18–23] is an ideal method due to the simplicity of its implementation as well as the direct explanations of various phenomena (including spatial-temporal effects, high-field transport, and ultrafast phenomena).

Therefore, we use an ensemble Monte Carlo method to study the far-field nonlinear THz radiation responses of the n-doped bulk InSb and InAs crystals, whose carrier concentrations are both about 10¹⁶/cm³ (See Fig. 1). It is found that due to the complex effects, the external radiation produces a great distortion in the time-domain spectrum, with new nonlinear frequency components used in the frequency spectrum. The three mechanisms, i.e., nonparabolicity, impact ionization and intervalley scattering, perform in different ways in different samples under the high THz field excitation. Partial results of the calculation are consistent with the experimental result of Ho and Zhang.^[17] The study has potential values in up-convension of THz wave and THz nonlinear frequency multiplication field.

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (a) Intense THz pulse induced far-field radiation in n-doped bulk samples at room temperature. (b) The illustration for the dynamic mechanisms. Electron energy deviates from that in parabolic band structure in the high THz field. Impact ionization process excites electron–hole pairs from the valence band. The intervalley scattering process causes the electrons to transfer from the Γ valley to satellite valleys.

We can use the Boltzmann equations to describe the dynamical evolution of the carriers in the THz field, which is studied by the ensemble Monte Carlo method as follows:^[20–23]

Carriers can have high kinetic energy E in high-field transport calculations, and the values of wave vector k are far from the minima of the conduction band. In this region, the energy deviates from the value calculated from the simple quadratic expression of a parabolic band structure, and nonparabolicity occurs. The E–k relationship is expressed as^[19]

The energetic electrons in the conduction bands (with high kinetic energy under strong electric field) will collide on lattice atoms and excite electron–hole pairs from the valence band. This phenomenon is called impact ionization. We apply the impact ionization mechanisms to our EMC model since they may play a very important role in narrow band gap semiconductors. There is a probability per unit time in which an impact ionization event happens, introduced by means of the Keldysh expression^[25]

We use 1 mm×1 mm material as a uniform system because the width of each direction is very large. Although this assumption will bring some errors, the small errors can be neglected. In the simulation, the conduction band structure is described in terms of the three-valley Γ–L–X model. The electric field, refined with a time step of 1 fs is assumed to be uniform in each time period. The number of electrons and the electric field are updated at each sampling time. The holes are supposed to be immobile, with electron transport taken into account. In addition, although the material volume is very large, we still consider the boundary. In the case of no bias, the Neumann boundary condition is used, namely when the particles reach the boundary, they will experience a mirror reflection process. Table 1 shows the parameters^[25–28] used for calculation.

Table 1.

Table 1.

Table 1. Parameters used for calculation. . Parameters InSb InAs Unit Density 5.790 5.667 g/cm3 Sound velocity 4060 4282 m/s Band gap 0.18 0.354 eV Lattice parameter 6.479 6.058 Å Intervalley Energy EΓ− L 0.76 0.73 eV EΓ− X 0.45 1.02 eV Effective mass mΓ 0.014 0.023 m0 mL(mX) 0.22(0.13) 0.29(0.64) m0 LO phonon energy 24.4 30 meV Intervalley phonon energy 19.9 17.4 eV Static dielectric function ɛ0 17.65 15.15 High-frequencydielectric function ɛ∞ 17.65 15.15

Table 1. Parameters used for calculation. .

The THz frequency is much lower than that of visible light, which is unable to directly excite semiconductor electrons from the valence band to the conduction band. However, free carriers directly accelerated by the THz field can be observed in n-doped semiconductors with high mobilities.^[29] For example, electrons can be accelerated into the ballistic regime in lightly doped GaAs, whose electron energies can be achieved to be larger than 1 eV in this way.^[29] A similar physical process is also easy to achieve for InSb and InAs, which have very high mobilities of ∼7.6×10⁴ cm²/(V·s) and ∼ 3.0×10⁴ cm²/(V·s), respectively. In many THz-pump/THz-probe experiments,^[15–17] quite a number of interesting experimental phenomena have been explained by using the principle of the THz field driven free electron. This theory has also been discussed by Dai et al.^[30] Therefore, we use the intense THz-field-driven electron to study the far-field radiation of the materials.

As is well known, the far-field terahertz radiation is actually determined by the current acceleration. In the strong driving field, the current acceleration is no longer linearly proportional to the driving field, making other components in the frequency domain occur in the outward radiation of THz wave.

Figure 2 shows the driving THz pulses of the frequency spectra and the time-domain spectra (the inset). Here the THz oscillation period is about 2.5 picosecond (ps), the center simulation frequency of the wave 1 THz, and the width 0.5 THz. According to the surface radiation calculations,^[31] the far-field terahertz formula is

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Driving THz pulses of the frequency spectra and the time–domain spectra (the inset). The THz oscillation period is about 2.5 ps, the center frequency of the wave 1 THz, and the width 0.5 THz.

Carriers are accelerated in the outfield, producing a variety of scatterings. The scattering rates and the different energy states show a complex relationship: they are related to the various scattering energy responses and their proportions.

Figure 3 shows the scattering rates versus energy in InSb (a) and InAs (b) at room temperature in the Γ valley. In the caculation, 13 scattering rates are considered, such as ionized impurity scattering, acoustic phonon scattering, intervalley scattering, intravalley scattering, polar optical phonon scattering and impact ionization. Nine relatively high scattering rates are demonstrated among them. The impact ionization and intervalley scattering become the major scattering mechanisms at high energy.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. Scattering rates versus energy in InSb (a) and InAs (b) at room temperature in the Γ valley. The different scattering rates are those of ionized impurity scattering (black line), acoustic phonon scattering (black dash curve), phonon absorption of Γ–L intravalley scattering (blue curve), phonon emission of Γ–L intravalley scattering (blue dash curve), phonon absorption of Γ –X intravalley scattering (green curve), phonon emission of Γ–X intravalley scattering (green dash curve), absorption (ab) of polar optical phonon scattering (pink curve), emission (em) of polar optical phonon scattering (pink dash curve), and impact ionization(red curve), respectively.

When a single intense THz pulse irradiates on the sample, free carriers are accelerated in the conduction band and the kinetic energy increases rapidly. In this region, the values of k are far from the minima of the conduction band, and the nonparabolicity occurs. Then, the rising strength of the THz field will lead to the impact ionization process because it requires lower energy than intervalley scattering for the narrow band gap semiconductors. The energetic electrons in the conduction bands, which have high kinetic energy, will excite electron–hole pairs from the valence band, thereby increasing the number of the carriers. When the intervalley scattering energy threshold is reached, the intervalley scattering process will take place.

Carriers dynamics has a direct influence on the far-field radiation, and the carrier distributions under different THz field strengths are calculated. We consider the maximum amplitudes of the incident THz field intensity in a range of 0 kV/cm–150 kV/cm. Figure 4 shows the electron percentage occupancies in different valleys and the addition percentage of electrons excited by impact ionization. Since the impact ionization generates new electron–hole pairs, the total number of carriers is updated at each time step. For InSb, impact ionization rises rapidly from 20 kV/cm, and the new particles (ΔN) increase up to about 87% of the total. The Γ–X intervalley scattering occurs from about 30 kV/cm, but no more than 5% of the final particles in the X valley can be reached, and almost no electrons are scattered into the L valley. Therefore, impact ionization is the dominant mechanism for n-doped InSb in the whole process. For InAs, the band gap (0.354 eV) is larger than InSb (0.18 eV), thus the impact ionization requires higher energy. Figure 4(b) shows that the impact ionization begins from 25 kV/cm, immediately after which the intervalley scattering occurs. Both impact ionization and intervalley scattering change slowly after the THz field has reached 100 kV/cm. As can be seen, for InAs, only about 38% of the adding particles can be reached in the impact ionization process, which is far below the percentage in the case of InSb, and 24% of the number of electrons in the L valley can be reached in the intervalley scattering process. The scattering mainly occurs in the Γ–L valley, and the intervalley scattering into the X valley is relatively less. Compared with Fig. 3(b), while the intervalley scattering has the higher ratio at the late stage, the impact ionization has the low threshold energy. The impact ionization will occur first with a higher proportion in the actual scatterings. As a result, both impact ionization and intervalley scattering work together in n-doped InAs.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Carriers distributions of InSb (a) and InAs (b) in different THz fields. The black, red, and blue curves are the percentages of electrons in the Γ, X, and L valleys respectively. The green curve is the percentage of new electrons excited by impact ionization.

This calculation result is consistent with the experimental result of Ho and Zhang.^[17] The work focused on the varying impurity doping types of bulk InSb and InAs materials. By using the THz pump-THz probe experiment, the reflection spectroscopy and nonlinear optical effects were studied in different THz fields (30, 63, 90, and 110 kV/cm). The results showed that for n-doped InSb, impact ionization dominated while for n-doped InAs, both impact ionization and intervalley scattering mechanisms had to be considered. In their discussion about the electron population, the changes of free carrier density at different THz field strengths, as well as the electrons transference to the X valleys were not considered. Studying the radiation of driving electrons under the intense THz fields, we obtain the same conclusion by the more accurate numerical simulation with more considerations.

We calculate the far-field radiations of the two materials according to formula (4) (See Fig. 5). The actual maximum amplitude of the incident THz field is chosen from 5 kV/cm–100 kV/cm. Equation (2) defines the concept of non-parabolic coefficient, directly determines the relationship between the coefficients of the electron energy and momentum. For InSb and InAs, since the band gap energy and the effective mass of an electron are relatively small compared with those of the other semiconductors, they have very high non-parabolic coefficients. Figures 5(b) and 5(d) show that at low field (peak 5 kV/cm), the shape of the radiation field is inconsistent with driven force (See Fig. 2). When the field is upraised, the nonlinear effect is subsequently generated. As a result, the radiation waveforms are distorted in the time–domain spectrum, and new frequency components are generated in the frequency spectrum. Since the non-parabolic coefficient of InSb is much higher than that of InAs, the two materials show different behaviors in certain THz outfields. With the peak amplitude on 20 kV/cm and 50 kV/cm, the new frequency peaks in InSb are located at about 3 THz and 5 THz, respectively. While for InAs, the radiation field is also in linear response to the driving field at 20 kV/cm, and the peak in the vicinity of 3 THz occurs till 50 kV/cm. In this drive strength, the new frequency peaks of the two materials are formed due to the band nonparabolicity. The intervalley scattering and impact ionization mechanisms of two materials will gradually appear, playing a leading role in the nonlinear effect when the outfield exceeds 50 kV/cm.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Far-field THz radiations under the excitations of different intense THz fields: 5 kV/cm (black curve), 20 kV/cm (red curve), 50 kV/cm (blue curve), and 100 kV/cm (green curve), respectively. Panels (a) and (b) show the frequency spectra and time-domain spectra for n-doped InSb; panels (c) and (d) display the frequency spectra and time-domain spectra for n-doped InAs.

In the high field of 100 kV/cm, the electron motion characteristics change significantly. Comparing the time-domain spectra in Figs. 5(b) and 5(d) with those in Fig. 2, the sharpening and asymmetry occur. For InSb, the peak increases dramatically after 1.5 ps, the impact ionization here plays a key role in addition to the nonparabolicity factor. On the one hand, the impact ionization can absorb a lot of electron kinetic energy, which slows the acceleration process; on the other hand, it can inspire additional carriers so that the current rise space enlarges. The effect of the latter is much greater than that of the former in the high field. For InAs in the field of 100 kV/cm, the peak of the radiation form is lower than those in the fields of 50 kV/cm and 20 kV/cm after 1.5 ps in Fig. 5(d). The main reason is that in addition to impact ionization, the intervalley scattering has a relatively high proportion, making the electrons transfer to the high energy valley. As a result, the effective mass increases rapidly while transforming a large initial kinetic energy. It can make the electronic group velocity plummet. Meanwhile, it can suppress the current upgraded by the impact ionization. Therefore, in the high THz field, the impact ionization and intervalley scattering make the process complex, leading to the great distortion of the InAs external radiation form.

The frequency spectrum shows that the radiation is further enhanced to about 7 THz under the high field of 100 kV/cm in InSb, whereas only 5 THz is generated in InAs. They both generate a 2-THz periodic regular distribution in the overall spectrum. The formations of these frequency peaks are closely related to the band nonparabolicity and the different scattering mechanisms of two materials. The reasons for the differences mainly lie in four aspects: firstly, the InSb nonparabolic coefficients are higher than InAs’s. Secondly, the effective mass of InAs is larger than InSb’s, and the electronic group velocity is relatively low. Thirdly, the impact ionization ratio of InSb is much higher than that of InAs. Thus more electrons are excited and make more contributions to high frequency. Fourthly, from 20 kV/cm to 100 kV/cm, the intervalley scattering ratio of InAs is higher than that of InSb, and the effective mass becomes large when electrons are scattered into the satellite valley, thus reducing the occurrences of high frequency components. The above factors all contribute to the new lower frequency of InAs than that of InSb.

In this work, we study the nonlinear radiation responses of two different n-type bulk semiconductors, i.e., InSb and InAs in a high THz field by using the ensemble Monte Carlo method. We also calculate the electron valley occupancy and the percentage of new electrons excited by impact ionization. The results show that firstly, the radiation frequencies in different-strength THz fields have different characteristics. Low field radiation is linear with respect to the outfield, while intense THz excitation beams trigger higher THz frequency radiations, indicating harmonic generation. The periodic regular distribution is present in the frequency spectrum. Moreover, there are different scattering mechanisms in different materials. The impact ionization dominates in n-doped InSb, while impact ionization and intervalley scattering work together in n-doped InAs. This result is consistent with the experimental result of Ho and Zhang.^[17] Finally, the peak value of each new frequency is actually due to the different nonlinear mechanisms. The band nonparabolicity and impact ionization promote the generation of nonlinear high frequency radiation, while intervalley scattering performs in the opposite way. The results of the work have potential values in up-convension of THz wave and THz nonlinear frequency multiplication field. In the simulation, the factors which do not influence the result analysis are normalized, and the results are relative values. Future calculations can be considered to improve the model to obtain the absolute value of the far-field THz radiation, so the energy conversion efficiency can be calculated. If the submicron-size materials are studied, the small-size effect and boundary conditions will result in the rapid change of the internal electric field and a series of much more complex nonlinear effects. In this case, the electric field should be updated by using the Poisson equation. In addition, we can also study the carrier dynamic effects of higher band gap materials like GaAs and CdTe in the high THz field by using this method.